package com.poxz.algorithm.template;

import java.util.*;

/**
 * 数学算法
 *
 * @author poxz
 */
public class MathOps {

    // ==================== 质数相关 ====================

    /**
     * 试除法判断质数
     * 时间复杂度: O(√n)
     *
     * @param n 待判断的数
     * @return 是否为质数
     */
    static boolean isPrime(int n) {
        if (n < 2) return false;
        for (int i = 2; i <= n / i; i++) {
            if (n % i == 0) return false;
        }
        return true;
    }

    /**
     * 埃氏筛法求质数（适合求1~n内所有质数）
     * 时间复杂度: O(n log log n)
     *
     * @param n 范围上限
     * @return 1~n内所有质数的列表
     */
    static List<Integer> sieveOfEratosthenes(int n) {
        boolean[] notPrime = new boolean[n + 1];
        List<Integer> primes = new ArrayList<>();
        
        for (int i = 2; i <= n; i++) {
            if (!notPrime[i]) {
                primes.add(i);
                // 筛掉i的所有倍数
                for (int j = i; j <= n / i; j++) {
                    notPrime[i * j] = true;
                }
            }
        }
        return primes;
    }

    /**
     * 欧拉筛（线性筛）- 最优质数筛法
     * 时间复杂度: O(n)
     * 每个合数只会被它的最小质因子筛掉一次
     *
     * @param n 范围上限
     * @return 1~n内所有质数的列表
     */
    static List<Integer> eulerSieve(int n) {
        boolean[] notPrime = new boolean[n + 1];
        List<Integer> primes = new ArrayList<>();
        
        for (int i = 2; i <= n; i++) {
            if (!notPrime[i]) {
                primes.add(i);
            }
            for (int prime : primes) {
                if (i * prime > n) break;
                notPrime[i * prime] = true;
                if (i % prime == 0) break; // 关键：保证每个数只被最小质因子筛掉
            }
        }
        return primes;
    }

    // ==================== 约数相关 ====================

    /**
     * 试除法求所有约数
     * 时间复杂度: O(√n)
     *
     * @param n 待求约数的数
     * @return 所有约数的有序列表
     */
    static List<Integer> getDivisors(int n) {
        List<Integer> res = new ArrayList<>();
        for (int i = 1; i <= n / i; i++) {
            if (n % i == 0) {
                res.add(i);
                if (i != n / i) {
                    res.add(n / i);
                }
            }
        }
        Collections.sort(res);
        return res;
    }

    /**
     * 求约数个数
     * 公式: 若 n = p1^a1 * p2^a2 * ... * pk^ak
     * 则约数个数 = (a1+1) * (a2+1) * ... * (ak+1)
     *
     * @param n 待求约数个数的数
     * @return 约数个数
     */
    static int countDivisors(int n) {
        int count = 0;
        for (int i = 1; i <= n / i; i++) {
            if (n % i == 0) {
                count++;
                if (i != n / i) count++;
            }
        }
        return count;
    }

    /**
     * 求约数之和
     * 公式: 若 n = p1^a1 * p2^a2 * ... * pk^ak
     * 则约数之和 = (1+p1+...+p1^a1) * (1+p2+...+p2^a2) * ... * (1+pk+...+pk^ak)
     *
     * @param n 待求约数之和的数
     * @return 约数之和
     */
    static long sumDivisors(int n) {
        long sum = 0;
        for (int i = 1; i <= n / i; i++) {
            if (n % i == 0) {
                sum += i;
                if (i != n / i) sum += n / i;
            }
        }
        return sum;
    }

    // ==================== 质因数分解 ====================

    /**
     * 试除法分解质因数
     * 时间复杂度: O(√n)
     *
     * @param n 待分解的数
     * @return 质因数及其次数的映射 Map<质因数, 次数>
     */
    static Map<Integer, Integer> primeFactorization(int n) {
        Map<Integer, Integer> factors = new LinkedHashMap<>();
        
        for (int i = 2; i <= n / i; i++) {
            if (n % i == 0) {
                int count = 0;
                while (n % i == 0) {
                    n /= i;
                    count++;
                }
                factors.put(i, count);
            }
        }
        if (n > 1) {
            factors.put(n, 1);
        }
        return factors;
    }

    // ==================== GCD 和 LCM ====================

    /**
     * 欧几里得算法求最大公约数（GCD）
     * 时间复杂度: O(log min(a,b))
     *
     * @param a 第一个数
     * @param b 第二个数
     * @return 最大公约数
     */
    static int gcd(int a, int b) {
        return b == 0 ? a : gcd(b, a % b);
    }

    /**
     * 最小公倍数（LCM）
     * 公式: lcm(a,b) = a * b / gcd(a,b)
     *
     * @param a 第一个数
     * @param b 第二个数
     * @return 最小公倍数
     */
    static long lcm(int a, int b) {
        return (long) a * b / gcd(a, b);
    }

    /**
     * 扩展欧几里得算法
     * 求解 ax + by = gcd(a,b) 的一组整数解
     *
     * @param a 第一个数
     * @param b 第二个数
     * @return 数组 [gcd(a,b), x, y]
     */
    static long[] exgcd(long a, long b) {
        if (b == 0) {
            return new long[]{a, 1, 0};
        }
        long[] res = exgcd(b, a % b);
        long gcd = res[0], x = res[2], y = res[1] - (a / b) * res[2];
        return new long[]{gcd, x, y};
    }

    // ==================== 快速幂 ====================

    /**
     * 快速幂 - 求 a^b mod m
     * 时间复杂度: O(log b)
     *
     * @param a 底数
     * @param b 指数
     * @param m 模数
     * @return a^b mod m
     */
    static long quickPow(long a, long b, long m) {
        long res = 1;
        a %= m;
        while (b > 0) {
            if ((b & 1) == 1) {
                res = res * a % m;
            }
            a = a * a % m;
            b >>= 1;
        }
        return res;
    }

    /**
     * 矩阵快速幂（用于求斐波那契数列等）
     * 时间复杂度: O(log n)
     *
     * @param n 第n项
     * @param m 模数
     * @return 斐波那契数列第n项 mod m
     */
    static long fibonacci(long n, long m) {
        if (n <= 0) return 0;
        if (n == 1 || n == 2) return 1;
        
        long[][] base = {{1, 1}, {1, 0}};
        long[][] res = matrixPow(base, n - 1, m);
        return res[0][0];
    }

    /**
     * 矩阵快速幂辅助函数
     */
    private static long[][] matrixPow(long[][] a, long n, long m) {
        int size = a.length;
        long[][] res = new long[size][size];
        // 初始化为单位矩阵
        for (int i = 0; i < size; i++) {
            res[i][i] = 1;
        }
        
        while (n > 0) {
            if ((n & 1) == 1) {
                res = matrixMultiply(res, a, m);
            }
            a = matrixMultiply(a, a, m);
            n >>= 1;
        }
        return res;
    }

    /**
     * 矩阵乘法
     */
    private static long[][] matrixMultiply(long[][] a, long[][] b, long m) {
        int n = a.length;
        long[][] c = new long[n][n];
        for (int i = 0; i < n; i++) {
            for (int j = 0; j < n; j++) {
                for (int k = 0; k < n; k++) {
                    c[i][j] = (c[i][j] + a[i][k] * b[k][j]) % m;
                }
            }
        }
        return c;
    }

    // ==================== 组合数 ====================

    /**
     * 递推法求组合数 C(n, m)
     * 杨辉三角: C(n,m) = C(n-1,m) + C(n-1,m-1)
     * 适用于多次查询，n ≤ 2000
     *
     * @param maxN 最大的n值
     * @param mod  模数
     * @return C[i][j] 表示 C(i,j) mod mod
     */
    static long[][] computeCombinations(int maxN, long mod) {
        long[][] c = new long[maxN + 1][maxN + 1];
        for (int i = 0; i <= maxN; i++) {
            c[i][0] = 1;
            for (int j = 1; j <= i; j++) {
                c[i][j] = (c[i - 1][j] + c[i - 1][j - 1]) % mod;
            }
        }
        return c;
    }

    /**
     * 通过定义式求组合数 C(n, m) = n! / (m! * (n-m)!)
     * 适用于 n ≤ 10000
     *
     * @param n   总数
     * @param m   选取数
     * @param mod 模数
     * @return C(n,m) mod mod
     */
    static long combination(int n, int m, long mod) {
        if (m > n || m < 0) return 0;
        if (m == 0 || m == n) return 1;
        
        // 预处理阶乘
        long[] fact = new long[n + 1];
        fact[0] = 1;
        for (int i = 1; i <= n; i++) {
            fact[i] = fact[i - 1] * i % mod;
        }
        
        // C(n,m) = n! / (m! * (n-m)!)
        // 在模意义下需要用逆元
        return fact[n] * quickPow(fact[m], mod - 2, mod) % mod 
               * quickPow(fact[n - m], mod - 2, mod) % mod;
    }

    /**
     * Lucas定理求大组合数 C(n, m) mod p （p为质数）
     * C(n,m) mod p = C(n/p, m/p) * C(n%p, m%p) mod p
     * 适用于 p ≤ 10^5, n,m ≤ 10^18
     *
     * @param n 总数
     * @param m 选取数
     * @param p 质数模数
     * @return C(n,m) mod p
     */
    static long lucas(long n, long m, long p) {
        if (m == 0) return 1;
        return combination((int) (n % p), (int) (m % p), p) 
               * lucas(n / p, m / p, p) % p;
    }

    // ==================== 其他常用算法 ====================

    /**
     * 求逆元（费马小定理）
     * 若p为质数，则 a^(p-1) ≡ 1 (mod p)
     * 因此 a^(-1) ≡ a^(p-2) (mod p)
     *
     * @param a 待求逆元的数
     * @param p 质数模数
     * @return a的逆元 mod p
     */
    static long modInverse(long a, long p) {
        return quickPow(a, p - 2, p);
    }

    /**
     * 欧拉函数 φ(n) - 求小于等于n且与n互质的数的个数
     * 公式: φ(n) = n * (1 - 1/p1) * (1 - 1/p2) * ... * (1 - 1/pk)
     * 其中 p1, p2, ..., pk 为n的所有质因数
     *
     * @param n 待求欧拉函数值的数
     * @return φ(n)
     */
    static int euler(int n) {
        int res = n;
        for (int i = 2; i <= n / i; i++) {
            if (n % i == 0) {
                res = res / i * (i - 1);
                while (n % i == 0) {
                    n /= i;
                }
            }
        }
        if (n > 1) {
            res = res / n * (n - 1);
        }
        return res;
    }

    /**
     * 欧拉函数筛法 - 线性求1~n的所有欧拉函数值
     * 时间复杂度: O(n)
     *
     * @param n 范围上限
     * @return phi[i] 表示 φ(i)
     */
    static int[] eulerSievePhi(int n) {
        int[] phi = new int[n + 1];
        boolean[] notPrime = new boolean[n + 1];
        List<Integer> primes = new ArrayList<>();
        
        phi[1] = 1;
        for (int i = 2; i <= n; i++) {
            if (!notPrime[i]) {
                primes.add(i);
                phi[i] = i - 1;
            }
            for (int prime : primes) {
                if (i * prime > n) break;
                notPrime[i * prime] = true;
                if (i % prime == 0) {
                    phi[i * prime] = phi[i] * prime;
                    break;
                } else {
                    phi[i * prime] = phi[i] * (prime - 1);
                }
            }
        }
        return phi;
    }

    /**
     * 判断完全平方数
     *
     * @param n 待判断的数
     * @return 是否为完全平方数
     */
    static boolean isPerfectSquare(long n) {
        if (n < 0) return false;
        long sqrt = (long) Math.sqrt(n);
        return sqrt * sqrt == n;
    }

    /**
     * 整数开方（向下取整）
     *
     * @param n 待开方的数
     * @return ⌊√n⌋
     */
    static long intSqrt(long n) {
        if (n < 0) return -1;
        long l = 0, r = n;
        while (l < r) {
            long mid = l + (r - l + 1) / 2;
            if (mid <= n / mid) {
                l = mid;
            } else {
                r = mid - 1;
            }
        }
        return l;
    }

    /**
     * 功能测试
     *
     * @param args 命令行参数
     */
    public static void main(String[] args) {
        System.out.println("=== 数学算法示例 ===\n");

        // 质数判断
        System.out.println("--- 质数判断 ---");
        System.out.println("17是质数? " + isPrime(17));
        System.out.println("18是质数? " + isPrime(18));
        System.out.println();

        // 质数筛
        System.out.println("--- 质数筛（1~30） ---");
        System.out.println("埃氏筛: " + sieveOfEratosthenes(30));
        System.out.println("欧拉筛: " + eulerSieve(30));
        System.out.println();

        // 约数
        System.out.println("--- 约数（n=12） ---");
        System.out.println("所有约数: " + getDivisors(12));
        System.out.println("约数个数: " + countDivisors(12));
        System.out.println("约数之和: " + sumDivisors(12));
        System.out.println();

        // 质因数分解
        System.out.println("--- 质因数分解（n=60） ---");
        Map<Integer, Integer> factors = primeFactorization(60);
        for (Map.Entry<Integer, Integer> e : factors.entrySet()) {
            System.out.println(e.getKey() + "^" + e.getValue());
        }
        System.out.println();

        // GCD 和 LCM
        System.out.println("--- GCD 和 LCM ---");
        System.out.println("gcd(12, 18) = " + gcd(12, 18));
        System.out.println("lcm(12, 18) = " + lcm(12, 18));
        System.out.println();

        // 快速幂
        System.out.println("--- 快速幂 ---");
        System.out.println("2^10 mod 1000000007 = " + quickPow(2, 10, 1000000007));
        System.out.println("3^100 mod 1000000007 = " + quickPow(3, 100, 1000000007));
        System.out.println();

        // 斐波那契数列
        System.out.println("--- 斐波那契数列 ---");
        long mod = 1000000007;
        System.out.println("fib(10) mod " + mod + " = " + fibonacci(10, mod));
        System.out.println("fib(50) mod " + mod + " = " + fibonacci(50, mod));
        System.out.println();

        // 组合数
        System.out.println("--- 组合数 ---");
        long[][] c = computeCombinations(10, 1000000007);
        System.out.println("C(5,2) = " + c[5][2]);
        System.out.println("C(10,3) = " + c[10][3]);
        System.out.println();

        // 欧拉函数
        System.out.println("--- 欧拉函数 ---");
        System.out.println("φ(12) = " + euler(12));
        System.out.println("φ(17) = " + euler(17));
        System.out.println();

        // 完全平方数
        System.out.println("--- 完全平方数 ---");
        System.out.println("16是完全平方数? " + isPerfectSquare(16));
        System.out.println("17是完全平方数? " + isPerfectSquare(17));
        System.out.println("⌊√100⌋ = " + intSqrt(100));
    }
}

